{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# M5 Forecasting Competition GluonTS Template\n",
    "\n",
    "This notebook can be used as a starting point for participating in the [M5 forecasting competition](https://www.kaggle.com/c/m5-forecasting-accuracy/overview) using GluonTS-based tooling.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standard imports\n",
    "\n",
    "First we import standard data manipulation libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import mxnet as mx\n",
    "from mxnet import gluon\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import json\n",
    "import os\n",
    "from tqdm import tqdm\n",
    "from pathlib import Path"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also define globally accessible variables, such as the prediction length and the input path for the M5 data. Note that `single_prediction_length` corresponds to the length of the validation/evaluation periods, while `submission_prediction_length` corresponds to the length of both these periods combined.\n",
    "\n",
    "By default the notebook is configured to run in submission mode (`submission` will be `True`), which means that we use all of the data for training and predict new values for a total length of `submission_prediction_length` for which we don't have ground truth values available (performance can be assessed by submitting prediction results to Kaggle). In contrast, setting `submission` to `False` will instead use the last `single_prediction_length`-many values of our training set as validation points (and hence these values will not be used for training), which enables us to validate our model's performance offline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "single_prediction_length = 28\n",
    "submission_prediction_length = single_prediction_length * 2\n",
    "m5_input_path=\"./m5-forecasting-accuracy\"\n",
    "submission=True\n",
    "\n",
    "if submission:\n",
    "    prediction_length = submission_prediction_length\n",
    "else:\n",
    "    prediction_length = single_prediction_length"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reading the M5 data into GluonTS\n",
    "\n",
    "First we need to convert the provided M5 data into a format that is readable by GluonTS. At this point we assume that the M5 data, which can be downloaded from Kaggle, is present under `m5_input_path`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calendar = pd.read_csv(f'{m5_input_path}/calendar.csv')\n",
    "sales_train_validation = pd.read_csv(f'{m5_input_path}/sales_train_validation.csv')\n",
    "sample_submission = pd.read_csv(f'{m5_input_path}/sample_submission.csv')\n",
    "sell_prices = pd.read_csv(f'{m5_input_path}/sell_prices.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start the data convertion process by building dynamic features (features that change over time, just like the target values). Here, we are mainly interested in the event indicators `event_type_1` and `event_type_2`. We will mostly drop dynamic time features as GluonTS will automatically add some of these as part of many models' transformation chains."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cal_features = calendar.drop(\n",
    "    ['date', 'wm_yr_wk', 'weekday', 'wday', 'month', 'year', 'event_name_1', 'event_name_2', 'd'], \n",
    "    axis=1\n",
    ")\n",
    "cal_features['event_type_1'] = cal_features['event_type_1'].apply(lambda x: 0 if str(x)==\"nan\" else 1)\n",
    "cal_features['event_type_2'] = cal_features['event_type_2'].apply(lambda x: 0 if str(x)==\"nan\" else 1)\n",
    "\n",
    "test_cal_features = cal_features.values.T\n",
    "if submission:\n",
    "    train_cal_features = test_cal_features[:,:-submission_prediction_length]\n",
    "else:\n",
    "    train_cal_features = test_cal_features[:,:-submission_prediction_length-single_prediction_length]\n",
    "    test_cal_features = test_cal_features[:,:-submission_prediction_length]\n",
    "\n",
    "test_cal_features_list = [test_cal_features] * len(sales_train_validation)\n",
    "train_cal_features_list = [train_cal_features] * len(sales_train_validation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then go on to build static features (features which are constant and series-specific). Here, we make use of all categorical features that are provided to us as part of the M5 data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "state_ids = sales_train_validation[\"state_id\"].astype('category').cat.codes.values\n",
    "state_ids_un , state_ids_counts = np.unique(state_ids, return_counts=True)\n",
    "\n",
    "store_ids = sales_train_validation[\"store_id\"].astype('category').cat.codes.values\n",
    "store_ids_un , store_ids_counts = np.unique(store_ids, return_counts=True)\n",
    "\n",
    "cat_ids = sales_train_validation[\"cat_id\"].astype('category').cat.codes.values\n",
    "cat_ids_un , cat_ids_counts = np.unique(cat_ids, return_counts=True)\n",
    "\n",
    "dept_ids = sales_train_validation[\"dept_id\"].astype('category').cat.codes.values\n",
    "dept_ids_un , dept_ids_counts = np.unique(dept_ids, return_counts=True)\n",
    "\n",
    "item_ids = sales_train_validation[\"item_id\"].astype('category').cat.codes.values\n",
    "item_ids_un , item_ids_counts = np.unique(item_ids, return_counts=True)\n",
    "\n",
    "stat_cat_list = [item_ids, dept_ids, cat_ids, store_ids, state_ids]\n",
    "\n",
    "stat_cat = np.concatenate(stat_cat_list)\n",
    "stat_cat = stat_cat.reshape(len(stat_cat_list), len(item_ids)).T\n",
    "\n",
    "stat_cat_cardinalities = [len(item_ids_un), len(dept_ids_un), len(cat_ids_un), len(store_ids_un), len(state_ids_un)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can build both the training and the testing set from target values and both static and dynamic features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gluonts.dataset.common import load_datasets, ListDataset\n",
    "from gluonts.dataset.field_names import FieldName\n",
    "\n",
    "train_df = sales_train_validation.drop([\"id\",\"item_id\",\"dept_id\",\"cat_id\",\"store_id\",\"state_id\"], axis=1)\n",
    "train_target_values = train_df.values\n",
    "\n",
    "if submission == True:\n",
    "    test_target_values = [np.append(ts, np.ones(submission_prediction_length) * np.nan) for ts in train_df.values]\n",
    "else:\n",
    "    test_target_values = train_target_values.copy()\n",
    "    train_target_values = [ts[:-single_prediction_length] for ts in train_df.values]\n",
    "\n",
    "m5_dates = [pd.Timestamp(\"2011-01-29\", freq='1D') for _ in range(len(sales_train_validation))]\n",
    "\n",
    "train_ds = ListDataset([\n",
    "    {\n",
    "        FieldName.TARGET: target,\n",
    "        FieldName.START: start,\n",
    "        FieldName.FEAT_DYNAMIC_REAL: fdr,\n",
    "        FieldName.FEAT_STATIC_CAT: fsc\n",
    "    }\n",
    "    for (target, start, fdr, fsc) in zip(train_target_values,\n",
    "                                         m5_dates,\n",
    "                                         train_cal_features_list,\n",
    "                                         stat_cat)\n",
    "], freq=\"D\")\n",
    "\n",
    "test_ds = ListDataset([\n",
    "    {\n",
    "        FieldName.TARGET: target,\n",
    "        FieldName.START: start,\n",
    "        FieldName.FEAT_DYNAMIC_REAL: fdr,\n",
    "        FieldName.FEAT_STATIC_CAT: fsc\n",
    "    }\n",
    "    for (target, start, fdr, fsc) in zip(test_target_values,\n",
    "                                         m5_dates,\n",
    "                                         test_cal_features_list,\n",
    "                                         stat_cat)\n",
    "], freq=\"D\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Just to be sure, we quickly verify that dataset format is correct and that our dataset does indeed contain the correct target values as well as dynamic and static features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "next(iter(train_ds))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the estimator\n",
    "\n",
    "Having obtained our training and testing data, we can now create a GluonTS estimator. In our example we will use the `DeepAREstimator`, an autoregressive RNN which was developed primarily for the purpose of time series forecasting. Note however that you can use a variety of different estimators. Also, since GluonTS is mainly target at probabilistic time series forecasting, lots of different output distributions can be specified. In the M5 case, we think that the `NegativeBinomialOutput` distribution best describes the output.\n",
    "\n",
    "For a full list of available estimators and possible initialization arguments see https://gluon-ts.mxnet.io/api/gluonts/gluonts.model.html.\n",
    "\n",
    "For a full list of available output distributions and possible initialization arguments see https://gluon-ts.mxnet.io/api/gluonts/gluonts.distribution.html."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from gluonts.mx import DeepAREstimator\n",
    "from gluonts.mx.distribution.neg_binomial import NegativeBinomialOutput\n",
    "from gluonts.mx.trainer import Trainer\n",
    "\n",
    "estimator = DeepAREstimator(\n",
    "    prediction_length=prediction_length,\n",
    "    freq=\"D\",\n",
    "    distr_output = NegativeBinomialOutput(),\n",
    "    use_feat_dynamic_real=True,\n",
    "    use_feat_static_cat=True,\n",
    "    cardinality=stat_cat_cardinalities,\n",
    "    trainer=Trainer(\n",
    "        learning_rate=1e-3,\n",
    "        epochs=100,\n",
    "        num_batches_per_epoch=50,\n",
    "        batch_size=32\n",
    "    )\n",
    ")\n",
    "\n",
    "predictor = estimator.train(train_ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generating forecasts\n",
    "\n",
    "Once the estimator is fully trained, we can generate predictions from it for the test values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from gluonts.evaluation.backtest import make_evaluation_predictions\n",
    "\n",
    "forecast_it, ts_it = make_evaluation_predictions(\n",
    "    dataset=test_ds,\n",
    "    predictor=predictor,\n",
    "    num_samples=100\n",
    ")\n",
    "\n",
    "print(\"Obtaining time series conditioning values ...\")\n",
    "tss = list(tqdm(ts_it, total=len(test_ds)))\n",
    "print(\"Obtaining time series predictions ...\")\n",
    "forecasts = list(tqdm(forecast_it, total=len(test_ds)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Local performance validation (if `submission` is `False`)\n",
    "\n",
    "Since we don't want to constantly submit our results to Kaggle, it is important to being able to evaluate performace on our own validation set offline. To do so, we create a custom evaluator which, in addition to GluonTS's standard performance metrics, also returns `MRMSSE` (corresponding to the mean RMSSE). Note that the official score for the M5 competition, the `WRMSSE`, is not yet computed. A future version of this notebook will replace the `MRMSSE` by the `WRMSSE`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "if submission == False:\n",
    "    \n",
    "    from gluonts.evaluation import Evaluator\n",
    "    \n",
    "    class M5Evaluator(Evaluator):\n",
    "        \n",
    "        def get_metrics_per_ts(self, time_series, forecast):\n",
    "            successive_diff = np.diff(time_series.values.reshape(len(time_series)))\n",
    "            successive_diff = successive_diff ** 2\n",
    "            successive_diff = successive_diff[:-prediction_length]\n",
    "            denom = np.mean(successive_diff)\n",
    "            pred_values = forecast.samples.mean(axis=0)\n",
    "            true_values = time_series.values.reshape(len(time_series))[-prediction_length:]\n",
    "            num = np.mean((pred_values - true_values)**2)\n",
    "            rmsse = num / denom\n",
    "            metrics = super().get_metrics_per_ts(time_series, forecast)\n",
    "            metrics[\"RMSSE\"] = rmsse\n",
    "            return metrics\n",
    "        \n",
    "        def get_aggregate_metrics(self, metric_per_ts):\n",
    "            wrmsse = metric_per_ts[\"RMSSE\"].mean()\n",
    "            agg_metric , _ = super().get_aggregate_metrics(metric_per_ts)\n",
    "            agg_metric[\"MRMSSE\"] = wrmsse\n",
    "            return agg_metric, metric_per_ts\n",
    "        \n",
    "    \n",
    "    evaluator = M5Evaluator(quantiles=[0.5, 0.67, 0.95, 0.99])\n",
    "    agg_metrics, item_metrics = evaluator(iter(tss), iter(forecasts), num_series=len(test_ds))\n",
    "    print(json.dumps(agg_metrics, indent=4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Converting forecasts back to M5 submission format (if `submission` is `True`)\n",
    "\n",
    "Since GluonTS estimators return a sample-based probabilistic forecasting predictor, we first need to reduce these results to a single prediction per time series. This can be done by computing the mean or median over the predicted sample paths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if submission == True:\n",
    "    forecasts_acc = np.zeros((len(forecasts), prediction_length))\n",
    "    for i in range(len(forecasts)):\n",
    "        forecasts_acc[i] = np.mean(forecasts[i].samples, axis=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then reshape the forecasts into the correct data shape for submission ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if submission == True:\n",
    "    forecasts_acc_sub = np.zeros((len(forecasts)*2, single_prediction_length))\n",
    "    forecasts_acc_sub[:len(forecasts)] = forecasts_acc[:,:single_prediction_length]\n",
    "    forecasts_acc_sub[len(forecasts):] = forecasts_acc[:,single_prediction_length:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ".. and verfiy that reshaping is consistent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if submission == True:\n",
    "    np.all(np.equal(forecasts_acc[0], np.append(forecasts_acc_sub[0], forecasts_acc_sub[30490])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we save our submission into a timestamped CSV file which can subsequently be uploaded to Kaggle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if submission == True:\n",
    "    import time\n",
    "\n",
    "    sample_submission = pd.read_csv(f'{m5_input_path}/sample_submission.csv')\n",
    "    sample_submission.iloc[:,1:] = forecasts_acc_sub\n",
    "\n",
    "    submission_id = 'submission_{}.csv'.format(int(time.time()))\n",
    "\n",
    "    sample_submission.to_csv(submission_id, index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting sample predictions\n",
    "\n",
    "Finally, we can also visualize our predictions for some of the time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plot_log_path = \"./plots/\"\n",
    "directory = os.path.dirname(plot_log_path)\n",
    "if not os.path.exists(directory):\n",
    "    os.makedirs(directory)\n",
    "    \n",
    "def plot_prob_forecasts(ts_entry, forecast_entry, path, sample_id, inline=True):\n",
    "    plot_length = 150\n",
    "    prediction_intervals = (50, 67, 95, 99)\n",
    "    legend = [\"observations\", \"median prediction\"] + [f\"{k}% prediction interval\" for k in prediction_intervals][::-1]\n",
    "\n",
    "    _, ax = plt.subplots(1, 1, figsize=(10, 7))\n",
    "    ts_entry[-plot_length:].plot(ax=ax)\n",
    "    forecast_entry.plot(prediction_intervals=prediction_intervals, color='g')\n",
    "    ax.axvline(ts_entry.index[-prediction_length], color='r')\n",
    "    plt.legend(legend, loc=\"upper left\")\n",
    "    if inline:\n",
    "        plt.show()\n",
    "        plt.clf()\n",
    "    else:\n",
    "        plt.savefig('{}forecast_{}.pdf'.format(path, sample_id))\n",
    "        plt.close()\n",
    "\n",
    "print(\"Plotting time series predictions ...\")\n",
    "for i in tqdm(range(5)):\n",
    "    ts_entry = tss[i]\n",
    "    forecast_entry = forecasts[i]\n",
    "    plot_prob_forecasts(ts_entry, forecast_entry, plot_log_path, i)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
